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ABSTRACT 
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Results  for  several  fighter  configurations  are  employed  for  comparison  with 
available  data.  It  is  shown  that  the  present  method  produces  reasonable 
results  in  predicting  sideslip  derivatives,  while  roll-  and  yaw-rate 
derivatives  do  not  compare  very  well  with  forced-oscillation  test  data  at  high 
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1.  INTRODUCTION 


Maneuvering  flight  at  high  angles  of  attack  may  induce  boundary-layer 
separation  and  vortex-separated  flow  over  aerodynamic  surfaces.  In  addition, 
the  usually  well-organized  vortex  flow  may  also  break  down  at  a  high  enough 
angle  of  attack.  These  flow  phenomena  affect  not  only  the  individual 
aerodynamic  surfaces  on  which  these  phenomena  occur  but  also  those  in  their 

close  proximity.  These  flow  patterns  are  now  well  known  qualitatively  for 

fighter  configurations  from  extensive  wind  and  water  tunnel  tests.  However, 
development  of  quantitative  prediction  methods  has  not  kept  pace  with  tunnel 
testing. 

Ultimately,  a  theoretical  method  for  predicting  aerodynamics  at  high 
angles  of  attack  would  be  ideally  based  on  Navier-Stokes  solutions  with 

appropriate  turbulence  modeling.  However,  for  applications  to  preliminary 
design,  this  is  not  feasible  at  the  present  time  due  to  inadequate  computer 
resources.  Therefore,  a  practical  approach  would  be  one  involving  panel  or 
panel-like  methods  with  corrections  for  high  angle-of-attack  flow  features. 
For  example,  one  currently  available  computer  code  for  complete  aircraft 
configurations  is  the  PAN  AIR  (refs.  1  and  2).  This  code  is  based  on  the 
solution  of  the  Prandtl-Glauert  equation  and  may  include  boundary-layer 
corrections.  However,  its  applications  to  configurations  with  extensive 
boundary-layer  separation  or  vortex  flow  have  not  been  demonstrated.  It 
should  be  noted  that  symmetrical  vortex  flow  without  breakdown  on  simple 
configurations  has  been  dealt  with  successfully  in  reference  3  based  on  a 

panel  method.  Similarly,  the  VSAERO  code  (refs.  4  and  5)  was  developed  mainly 
for  symmetrical  flight  conditions  with  vortex  flow  and  boundary-layer 
correction  at  moderate  angles  of  attack.  A  comprehensive  method  for  lateral- 
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directional  aerodynamics  with  aeroelastic  effect  is  that  of  the  FLEXSTAB  (ref. 

6) .  However,  it  was  valid  only  in  the  region  of  linear  aerodynamics,  and 
hence  at  low  angles  of  attack.  Recently,  the  VSAERO  code  has  been  extended  to 
enable  calculation  of  stability  and  control  characteristics  of  airplanes  (ref. 

7) .  Again,  applications  to  high  angle-of-attack  aerodynamics  have  not  been 
reported. 

To  remedy  the  inadequacy  of  existing  methods,  particularly  in  the  area  of 
nonlinear  lateral-directional  aerodynamics,  the  VORSTAB  code  was  developed  for 
vortex-dominated  configurations  (refs.  8  and  9).  In  this  code,  the  effect  of 
vortex  lift  is  included  through  the  method  of  suction  analogy  while  the  effect 
of  vortex  breakdown  has  been  accounted  for  using  a  theoretical  correlation 
parameter  obtained  from  experimental  data.  The  code  can  treat  conventional 
fighter  configurations  by  using  nonlinear  section  data  for  the  effect  of 
boundary-layer  separation  (ref.  10).  Either  could  be  used  in  the  code 
experimental  or  theoretical  nonlinear  section  data  in  order  to  avoid  extreme 
difficulty  in  the  theoretical  prediction  of  three-dimensional  viscous 
separated  flows. 

In  the  present  investigation,  the  method  of  suction  analogy  used  in 
reference  8  is  replaced  with  the  method  of  free  vortex  filaments  (ref.  11)  to 
allow  a  more  accurate  calculation  of  interference  effects  among  different 
lifting  surfaces.  In  addition,  forebody  vortex  separation  on  noncircular 
cross  sections  is  also  modeled.  The  results  obtained  demonstrate  the  promise 
of  the  present  code  on  a  design  tool. 

In  this  report  both  the  theoretical  and  empirical  methods  used  in  the  code 
will  be  described.  Numerical  results  to  demonstrate  the  applicability  of  the 
present  code  to  modern  fighter  aircraft  geometries  will  be  presented. 
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2.  THEORETICAL  APPROACHES 


2.1.  General  Concept 

The  main  flow  features  accounted  for  in  the  present  method  are 
illustrated  in  Figure  1.  The  basic  flow  considered  in  the  code  is  the 
subsonic  potential  flow,  with  corrections  for  boundary-layer  separation  and 
imbedded  free  vortex  filaments  for  vortex  flow  effect.  The  potential  flow 
calculation  is  based  on  the  small  disturbance  approximation  of  the  gas  dynamic 
equations  which  result  in  the  Prandtl-Glauert  equation.  There  are  two  kxuds 
of  flow  singularities  used  in  the  program:  horseshoe  vortices  for  lifting 
surfaces  and  vortex  multiplets  for  the  fuselage.  The  strength  of  the 
singularities  can  be  obtained  by  satisfying  the  boundary  conditions  on  the 
fuselage  and  lifting  surfaces.  To  properly  account  for  the  leading-edge 
singularity  of  pressure  loading  in  the  linear  theory,  and  hence  the  leading- 
edge  thrust,  the  Quasi-Vortex-Lattice  Method  (QVLM)  is  used  (ref.  12). 

In  addition,  the  vortex  flow  arising  from  separation  along  the  leading 
and  side  edges  of  a  low  aspect  ratio  wing  which  significantly  affects  the  high 
angle-of-attack  aerodynamics  is  simulated  by  using  discrete  free  vortex 
filaments  emanating  from  the  edges.  At  a  high  enough  angle  of  attack,  the 
vortex  flow  will  breakdown  and  decrease  the  vortex  lift  effect.  The  latter  is 
calculated  in  the  present  method  by  using  empirical  formulas  derived  from  the 
analysis  of  experimental  data. 

Another  phenomenon  limiting  the  high  angle-of-attack  flight  envelopes 
for  high  performance  aircraft  is  the  boundary-layer  separation.  To  account 
for  this  viscous  effect  on  the  lifting  surfaces,  sectional  nonlinear  data  are 
used  iteratively.  In  the  following,  the  various  methodologies  used  in  the 
code  will  be  summarized. 
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2.2.  Potential  Flow  Theory 


As  indicated  earlier,  the  present  potential  flow  method  is  based  on  the 
solution  of  the  Prandtl-Glauert  equation: 


2  2  2 

n  m2.  a  $  ,  o  ,  o 

dx  dy  dz 


(1) 


where  $  is  the  perturbation  velocity  potential.  The  solution  is  represented 
by  vortex  distributions.  The  thin  wing  approximation  is  used  throughout. 


2.2.1.  Boundary  Conditions 

Lifting  Surfaces.  The  boundary  condition  for  the  Prandtl-Glauert 
equation  on  the  wing  is  that  the  velocity  component  normal  to  the  wing  surface 
should  be  zero.  Assuming  that  the  wing  surface  can  be  described  as 

z  -  zc(x,y)  (2) 


The  unit  normal  vector  on  the  wing  surface  can  be  written  as 

j  +  £ 


dz  dz 

_ C  t _ C 

dx _ 9y 


/S'  2  _  « 

i  +  (^)2  +  (^)2 


(3) 


dz 
_ < 

■dy 


For  an  airplane  having  angular  motions  with  rates  p,  q,  and  r  being 
positive  along  the  positive  x',  y'  and  z'-direction,  respectively  (see  fig.  2 
for  coordinate  systems),  the  onset  velocity  at  (x',y',z')  relative  to  the 
airplane  due  to  these  angular  rates  is 
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“(ip+jq+kr)  x  (ix'+jy'+kz')  =*  -i(qz'-ry')  +  j(pz'-rx')  -  k(py'-qx') 


To  express  this  velocity  vector  in  the  aerodynamic  coordinate  system  (x,y,z), 
the  sign  of  i  and  k-  components  must  be  reversed  and  x',  y'  and  z'  are  to 
be  replaced  by  -x,  y,  and  -z,  respectively.  It  follows  that 

“►  •>  ■> 
v  =  i(-qz  -  ry)  +  j(-pz  +  rx)  +  k(py  +  qx) 

This  must  be  added  to  the  induced  velocity  (u,v,w)  due  to  a  vortex 
distribution  and  the  uniform  freestream  contribution  due  to  angle  of  attack 
(a)  and  sideslip  (p)  to  obtain  the  total  velocity  vector  at  (x,y,z): 

V  =  (V^  cos  a  -  qz  -  ry  +  u)i  +  (-V^  sin  p  -  pz  +  rx  +  v)j 

+  (Vro  sin  a  +  py  +  qx  +  w)k  (4) 

Applying  the  tangency  condition 

V  •  n  -  0  (5) 

the  boundary  condition  becomes 

5z  3z 

w  " v  *  (v» cos  a  _  qz  -  ry) 

dz 

+  (-V^  sin  p  -  pz  +  rx)  sin  a  -  py  -  qx  (6) 

where  udz  /5x  has  been  ignored  as  a  second-order  term.  For  a  wing  with 
c 

dihedral  ( 4> )  and  twist  (a  ) ,  it  can  be  shown  that  (ref.  13) 

tw 
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-  sin  a...  +  -z —  cos  a 

tw  dx2  tw 


<))(cos  a  +  ~z  *  sin  a  1 
YV-  tw  dx„  tw; 


sin  $  cos  a _  +  —  cos  <b 
tw  dy  2  Y 

dz 

cos  4>(cos  atw  -  ~  sin  4.) 


where  the  (X2,y2)  coordinate  system  is  illustrated  in  Figure  3.  Basically,  it 
is  a  coordinate  system  based  on  the  local  wing  chord  plane  in  a  nonplanar 
configuration.  For  a  cambered  wing  with  small  spanwise  slopes  (dzc/dy2)>  it 
can  be  shown  from  equation  (7b)  that 


Therefore,  the  wing  boundary  condition  becomes 


w  v  c 

V  =  —  cos  4»  -  —  sin  4>  ■  V  +  V  c  =  cos  4>  -r —  cos  a  -  sin  a  cos 
n  V  V  nw  nf  ox 


-  sin  p  sin  <(>  -  (-§~- )  ^  (z  sin  $  +  y  cos  <j>) 


'2V  1  b 

CO 


+  (fv”)  |  (*  sin  ♦  -  y  cos  4> 


2V  >  b 


—  0  oz 

~  (fv-)  =  (x  +  2  cos  0 


where  Vnw  is  the  induced  normal  velocity  due  to  the  vortex  distribution  of 
lifting  surfaces  and  is  that  due  to  the  fuselage  and  c  in  the  mean 
aerodynamic  chord. 
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Fuselage . 


The  fuselage  surface  will  be  described  with  the  following 


expression: 


F(r,0,x)  *  r  -  R(0,x) 


(10) 


where  the  coordinate  system  is  shown  in  Figure  4.  It  follows  that  the  normal 

vector  on  the  surface  in  cylindrical  coordinate  system  with  e  ,  e  as  basis 

r  0 

vectors  is  given  by 


VF 


-*  l  3R _ +  dR 

V  R  30'  "  1  3x 


1  dR 

er  +  e0'  R  50 


t 

dx 


(ID 


The  total  velocity  vector  is  still  given  by  equation  (4).  As  shown  in  Figure 

3  the  unit  vectors,  j  and  k,  are  related  to  the  basis  vectors,  e  and  e_, 

J  r  0 

in  the  following  way: 


3 

->■ 

■ 

sin  0 

-cos  0 

-y 

e 

r 

-► 

k 

= 

cos  0 

sin  0 

_v 

(12) 


Therefore,  the  total  velocity  vector  can  be  rewritten  as 


V  »  (V  cos  a  -  qz  -  ry  +  u)i  +  [v  sin  0  +  w  cos  0  +  rx  sin  0  +  qx  cos  e 

+  p(y  cos  0  -  z  sin  0)  -  Va  sin  p  sin  0  +  sin  a  cos  0]er 

+  [v  cos  0  -  w  sin  0  +  rx  cos  0  +  qx  sin  0  +  p(y  sin  0  +  z  cos  0) 

+  V  sin  p  cos  0  +  V  sin  a  sin  0]e_,  (13) 

00  00  J  0 
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Applying  the  boundary  condition: 


V  .  7F  =  0, 


it  is  obtained  that 


9  1  9  R  9  R  ,  ,  ,  1  9R 

v”  R  ae*  ‘  cos  “  sx  "  sln  “  cos  9  •  sln 


1  9  R 

+  sin  g  sin  0  -  sin  g  cos  9  —  — 


,  pb  r  y  cos  9  -  z  sin  9  y  sin  9  +  z  cos  9  1_  9R-. 
V2V  J  L  b/2  +  h/2  R 


rb__,  r  x  sin  9  x  cos  9  1.  9R.  y  9R-. 


+  r  -i  r  .  oj-“ 

^2V  b/2 


b/2  R  90  b/2  9xJ 


where 


,  qc>  rx  cos  9  x  sin  9  1  9R  z  9R, 

lzVL  dl  c/2  R  90  ”  c/2  3x 


vr  =  v  sin  9  +  w  cos  0 


=  v  cos  9  -  w  sin  0 


It  should  be  noted  that  the  term  associated  with  (pb/2V  )  can  be  simplified 

QO 

because  with  y  ■  R  sin  0  and  z  ■  R  cos  9 , 


y  cos  9  -  z  sin  0=0 


(15a) 


y  sin  9  +  z  cos  0  =  R 


(15b) 
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•  Vortex  Distributions  on  Lifting  Surfaces 


To  satisfy  the  Prandtl-Glauert  equation  (1)  and  the  boundary  condition 
(eq.  9),  surface  vortex  distributions  on  lifting  surfaces  and  the  wake  are 
used  to  produce  an  induced  velocity  vector  at  any  point  (R)  in  the  field 
given  by 


(16) 


(17) 


and  u)  is  the  vorticity  vector  at  =  x^i  +  y^j  +  z^k.  Equation  (16)  is  to 
be  substituted  into  equation  (9)  to  obtain  an  equation  for  the  unknown  u.  To 
simplify  the  arithmetic,  a  typical  procedure  is  to  reduce  the  integral  in 
equation  (16)  to  a  finite  sum  through  a  panel  method  or  a  vortex-lattice 
method.  In  the  present  code,  a  method  called  the  quasi-vortex-lattice  method 
(QVLM)  is  used  to  sum  the  effect  of  horseshoe  vortices  (ref.  12).  The  QVLM 
methodology  was  developed  to  accurately  account  for  mathematical  singularities 
of  the  square-root  type  at  the  leading  edges  and  the  Cauchy  type  in  the  chord- 
wise  integral.  Accurate  calculation  of  the  leading-edge  thrust  is  needed  in 
calculating  the  vortex-flow  effect  (see  section  2.3.1). 

To  illustrate  the  essential  idea  of  the  QVIii,  consider  a  thin  airfoil  in 
a  two-dimensional  flow.  The  induced  velocity  at  x  on  the  chord  line  is 

w(x,0)  -  -  ^  (18) 
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The  integrand  in  equation  (18)  possesses  a  Cauchy  singularity  at  £  *  x.  In 
addition,  y(£)  has  a  square-root  singularity  at  £  -  0.  To  eliminate  the 
square-root  singularity,  the  x-coordinate  is  transformed  to  a  0-coordinate 
through  the  following  relation: 


(1  -  cos  9) 
2 

(1  -  cos  9") 
2 


(19) 

(20) 


Therefore,  equation  (18)  becomes 


Let 


w(9) 


1_  rU  y(8")sin  9'  d 9 
2 it  ^  cos  0  -  cos  0' 


g(9)  =  y(0)  sin(9) 


(21) 

(22) 


Since  y(9)  has  a  square-root  singularity  at  9-0  and  sin  9  vanishes  at 
0*0  as  fx,  the  square-root  singularity  of  y(9)  is  eliminated  by  the 
factor  sin  0.  Equation  (21)  can  be  written  as 


v(9) 


JL  r%  g(9~)  -  g( 9)d9 
2x  ^  cos  9  -  cos  Q' 


(23) 


Equation  (23)  can  now  be  reduced  to  a  finite  sum  through  the  midpoint 
trapezoidal  rule  as  follows : 
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x^  is  x  at  0^  -  in/N 
x^  is  at  8^  -  (2k  -  l)x/2N 

where  C  is  related  to  the  frequently-used  leading-edge  singularity  parameter 
C' 

C'  -  lim  u(x)  /x  -  lim  -  y(x)/x  -  |  (25) 

x+0  x+0  2  1 

and  can  be  computed  once  y^'s  are  obtained. 

In  a  three-dimensional  flow,  the  wing  surface  with  a  leading  edge  sweep 
A^  is  divided  into  vortex  strips.  Over  each  vortex  strip,  the  vortex 
integral  is  reduced  to  a  finite  sum  in  exactly  the  same  manner  as  in  the  thin 
airfoil  case  described  above  with  the  leading-edge  singularity  parameter  C 
calculated  from  the  following  relation  (ref.  12): 

NC/tan^  +  (3^  *  \  (induced  upwash  at  Jl.e.) 
dzc 

-  [cos  <p  cos  a  -  cos  (J>  sin  a]  (26) 


The  sectional  leading-edge  thrust  coefficient  (cfc)  is  then  given  by 


iiC2/  1  -  M2  cos A, 


2  cos  A„ 


Detailed  expressions  for  the  Induced  velocity  vector  to  be  used  in  the 
boundary  conditions  can  be  found  in  reference  12. 
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2.2.3.  Distribution  of  Vortex  Multiplets  for  the  Fuselage  Effect 


For  the  fuselage  effect,  G.  N.  Ward's  vortex  multiplets  are 
distributed  along  the  fuselage  axis  (ref.  14).  The  velocity  potential  is 
given  by 


y-  *f  (x»9  ,r)  = 


^1__  Y  r  cos  n0  i  . 
4ir  4  t  sin  n0  ‘  J 


ft 


n 


ft 


.u-^x-o2+-j2Eir  f  (c)d5 

r11/  (x-^)2  +  62r2 


(28) 


where  cos  n0  is  for  the  symmetrical  loading  and  sin  ne  for  the 

antisymmetrical  loading.  In  an  asymmetrical  flow,  such  as  in  a  lateral  or 

directional  motion,  both  terms  are  needed.  To  use  equation  (28),  the  integral 

is  first  integrated  by  parts  and  the  resulting  expression  for  is  then 

differentiated  to  obtain  V  and  V  as  follows: 

r  0 


v  ,  3$.  ,  3F 

r  1 _  _ f  1  r  r  cos  n0 1  <  ,  .  .  _ n  .  . 

U  )f  V  3r  4tt  l  1  sin  neHfn(xft>  ar  (*ft> 


n 


3F 

f  (x_  )  — —  (x  ) 
n  f£  3r  v  f£ 


ft  3 F 


+  f 


(29) 


ft 


(hf 


V  ra© 


-n  sin  n0 

jt?h 

n  n  cos  n0 


f  (x.  )F  (xc  ) 
nv  fz  nv  fz 


ft 


+  /  Fn<6)f;(e>*} 


(30) 


ft 


where 
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n  *  0 


When  equations  (29)  and  (30)  are  substituted  into  the  boundary  condition  (eq. 
(14)),  the  resulting  equation  becomes  an  integral  equation  for  the  unknown 
singularity  strength  f'(£).  The  integrals  involving  f'(£)  are  reduced  to 
finite  sums  as  follows: 


ft  SF 

9r~  fn  d^  -  (xft  ~  xfJp  2NL  J  fnk  sin  9l 
xfl  F  k 


ft 

XfJL 


(xft  Xfl)  2N  ^  fnk  Sin  9k  Fnk 
F  k 


(33) 


(34) 


The  transformation  used  between  5  and  0  is  to  transform  £  »  (xf^,  xft) 
into  0  ■  (0,  ti/2). 

The  boundary  condition  (eq.  14)  is  to  be  satisfied  on  the  surface.  For 
this  purpose,  the  method  of  weighted  residuals  will  be  used  with  R^^cos  m0 
as  the  weighting  function,  where  R^(x^,0)  is  the  radius  function  at  a 
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control  station 


x^.  Now  equation  (14)  is  multiplied  by  cos  m0  and 

integrated  from  —it  to  it,  where  m=0,l,»«»  .  For  example,  in  a  symmetrical 
case  with  m*0: 


.  it  v  it  v  SR.  it  SR 

~  l  rp-  R,d0  -  it  I  T-  R^de  "  cos  o  /  R, 

2it  J  V  i  2 it  1  V  R.  S0  i  2it  1  i  Sx 

—it  00  —it  00  i  —it 


d0 


-  sm 


1  it  SRi  it 

a  J  sin  0  -r—  —  R  d0  -  sin  a  ■=—  /  cos  0  R.  d0, 

2it  1  R.  S©  i  2 it  J  l 

—it  i  —it 


!=!,•••, Np.  (35) 


Similarly,  for  m=l. 


,  it  v  it  v  SR 

u  I  T-  Ei  “s  Me  -  25  /  V-  S7  5T  Ri  “s  M9 

-it  oo  —it  00  i 


cos 


a 


-it 


R 


2  SR 
i  Sx 


cos  0d0  -  a  j  cos  0 

7ir  J 


It 


2lt 


-It 


cos  0d0 


sin  a  r 

2n  •* 
—it 


SR 


sin  0 


i  „2 


R.  30 


cos  0d0, 


(36) 


Note  that  v  and  v  include  both  fuselage  and  wing  effects.  Equations 

(35),  (36),  etc.  are  used  to  solve  for  the  unknowns  f'.  For  a  fuselage  with 

n 

finite  slope  at  the  nose,  f  (x, „)  =  0. 

n  fJL 

It  should  be  noted  that  when  a  combination  of  surface  and  axial  vortex 
distributions  is  used,  the  wing  loading  may  become  too  low  at  the  wing- 
fuselage  juncture  due  to  the  horseshoe-vortex  discretization  of  wing  vortex 
distributions.  This  problem  is  solved  in  the  present  code,  under  the 
symmetrical  loading  condition,  by  eliminating  the  inboard  trailing  vortices 
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associated  with  the  root  vortex  strip.  If  this  is  not  done,  the  root  vortex 
strip  will  exhibit  a  loading  similar  to  that  near  the  wing  tip. 

A  more  detailed  derivation  of  equations  for  a  body  of  revolution  can  be 
found  in  reference  15. 

2.3.  Vortex  Flow  Theory 

When  a  lifting  surface  is  assumed  to  exhibit  edge -separated  vortex  flow, 
free  vortex  filaments  emanating  from  these  edges  are  employed  to  simulate  the 
induced  aerodynamic  effect  (ref.  11).  The  strength  of  these  free  vortex 

filaments  is  affected  by  the  leading-edge  radius,  wing  thickness,  boundary 
layer  separation,  and  vortex  breakdown.  In  addition,  forebody  vortices  may 
produce  significant  effect  on  overall  aerodynamic  characteristics  of  a 

configuration.  Theoretical  methods  to  calculate  these  effects  are  discussed 
in  the  following. 

2.3.1.  Free  Vortex  Filaments  (Reference  11) 

The  idea  of  using  free  vortex  filaments  for  leading  edge  separated 

flow  was  described  in  reference  11.  A  typical  vortex  filament  is  shown  in 
Figure  5.  A  filament  consists  of  a  series  of  short  straight  segments.  The 
initial  location  of  these  segments  is  shown  in  dashed  lines  and  the  final 

location  in  solid  lines.  Points  A  through  E  lie  along  a  wing  trailing  vortex 
element.  BC  is  allowed  to  move  only  in  the  free  stream  direction,  whereas  CD 
is  fixed  in  the  wing  plane  because  the  flow  is  tangential  to  the  trailing 
edge.  DE  is  also  fixed  in  the  wing  plane.  Both  segments  AB  and  IJ  are 
allowed  to  move  to  align  with  the  local  velocity  vector. 
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In  addition  to  the  force-free  condition  on  all  free  vortex  elements, 
boundary  conditions  *re  also  satisfied  along  the  leading  edges.  These 
conditions  are  related  to  the  degree  of  vortex  separation.  Note  that  the 
amount  of  residual  leading-edge  thrust  can  be  used  as  a  measure  of  the  degree 
of  vortex  separation.  In  other  words,  if  the  leading-edge  thrust  is  zero, 
then  a  complete  separation  has  occurred.  Otherwise,  only  partial  separation 
exists.  Now,  it  is  known  that  the  leading-edge  thrust  is  related  to  the  net 
upwash  at  the  leading  edge  (ref.  12)  so  that  the  leading-edge  boundary 
condition  can  be  written  as 


V  =  a  +  V 
n  Z  «n 


(37) 


where  V  ^  is  given  by  the  right-hand  side  of  equation  (9)  and  a^  is  given 
by  (refs.  11  and  12) 


a  =  total  induced  normalwash  -  V 
Z  oo  n 


- —  2  Ct  cos  A  1/2 

=  N  /tan  A  +  6  ( — — - . ) 

Z  7  2  2 

tt/1  -  M  cos  A. 


(38) 


For  a  case  with  complete  leading-edge  separation,  c*.  ■  0  and  so  is  a  .  Note 

c  Z 

that  ct  may  be  zero  at  a  portion  of  the  leading  edge  and  not  zero  at  other 
parts,  depending  on  such  factors  as  leading-edge  radius  and  thickness 
distributions. 

Equation  (37),  together  with  the  force-free  conditions  on  all  free  vortex 
segments,  is  solved  iteratively  (ref.  11).  To  reduce  the  computing  time,  the 
total  number  of  vortex  segments  used  may  be  reduced  by  half  during  the  first 
few  iterations.  In  addition,  to  avoid  numerical  instability  associated  with 
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1 


crowded  vortex  filaments,  an  option  is  provided  to  calculate  the  induced 
velocity  at  midpoints  between  two  adjacent  filaments  and  interpolate  the 
results  for  the  induced  velocity  at  the  filaments  for  the  purpose  of 
satisfying  the  force-free  condition. 

2.3.2.  Effect  of  Leading-Edge  Radius  on  Vortex  Separation 

As  indicated  earlier,  the  strength  of  vortex  filaments  is  affected  by 
the  leading-edge  radius.  A  simple,  yet  effective  concept  is  needed  to  connect 
the  two.  In  the  present  method,  this  is  provided  by  Kulfan's  assumptions 
(refs.  16  and  17)). 

Kulfan  assumed  that  on  a  slender  wing  the  It iding-edge  vortex 

separation  starts  at  an  angle  of  attack  (as)  at  w^ich  the  leading-edge  drag 

equals  the  leading-edge  thrust  (refs.  16  and  17).  This  condition  has  been 

shown  to  be  consistent  with  the  angle  of  attack  at  which  the  leading-edge 

laminar  separation  first  occurs  (ref.  18).  To  calculate  a  ,  consider  a 

s 

cambered  wing.  The  sectional  leading-edge  thrust  coefficient  can  be  written 
as 

ct  =  K(sin  a  +  aQ (39) 

where  K  is  a  function  of  geometry  and  Mach  number,  and  is  the  sectional 

angle  of  zero  lift.  For  a  wing,  the  sectional  leading-edge  thrust  coefficient 
was  given  in  equation  27.  The  leading-edge  singularity  parameter  at  any  a 
can  be  written  as 

C  *  K'(sin  a  +  a  )  (40) 

o 

Therefore,  the  leading-edge  singularity  parameter  at  ag  is  obtained  as 
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(41) 


C.  -  C(sin  a  +  a  )/(sin  a  +  a  ) 
1  so  o 


The  starting  angle  of  attack  (ag)  of  the  leading-edge  vortex  separation  can 
be  obtained  by  equating  the  leading-edge  drag  to  the  leading-edge  thrust. 
Using  the  expression  for  the  leading-edge  drag  from  reference  19,  it  is 
obtained  that 


cos  A, 


2  22  1/2 
0  (sin  a  +  a  )  (1  -  M  cos  A„) 

±  m  71  r2  s  0  00  1  (A?} 

2  2 ,  vl/2  2  L  .  ,  a  cos  A„  ^  ' 


(1  -  M  cos  A0) 


(sin  a  +  a  ) 
o 


where  rQ  is  the  leading-edge  radius.  It  follows  that 


a  »  sin  s 
s  L 


With  a  calculated,  the  sectional  thrust  coefficient  at  a  >  a  is  given  by 
s  s 


c  =  (x/2)C^( 1  -  M2cos2A.)1/2/cos  A. 
t  i.  ®  x.  x. 


(44) 


where 


C9  =  C[sin(a  -  ag)  +  aQ]/(sin  a  +  an) 


(45) 


This  sectional  thrust  coefficient  is  subtracted  from  the  full-thrust  value  to 
obtain  the  residual  thrust  which  is  used  in  equation  (38). 


2.3.3.  Effect  of  Wing  Thickness  Distribution 

Experiments  show  that  the  leading-edge  vortex  cores  move  upwards  and 
outwards  due  to  thickness  (ref.  20).  This  may  be  due  to  the  following 
reasons.  Firstly,  wing  thickness  affects  the  induced  velocity  in  the  flow 


-18- 


field.  Secondly,  following  Kulfan  (ref.  16),  the  flow  from  the  lower  surface 
has  a  smaller  angle  to  turn  around  the  leading  edge  to  the  upper  surface  of  a 
thick  wing.  Therefore,  the  net  centrifugal  force  necessary  to  turn  the 
attached  flow  is  less  than  that  required  for  a  thin  wing.  It  follows  that  the 
leading-edge  thrust  is  reduced. 

To  account  for  the  induced  flow  produced  by  wing  thickness,  a  source 
distribution  is  used  with  the  velocity  potential  given  by 

_ g(x~,  y^,z~)ds _  (46) 

/(x-x')2  +  s2(y-y')2  +  e2(z-z')2 

where  a  is  the  source  strength: 


$  =  - 


_1 

4-n 


II 

S 


w 


a 


=  _  ril-j 

'•3n-,+  l-3nJ- 


(47) 


With  the  unit  normal  vector  due  to  thickness  given  by 


and 


* 


cos  a  1  +  sin  cx£. 


(48) 


the  outflow  velocity, 


i! 

3  n* 


can  be  shown  to  be 


(49) 


This  outflow  velocity  is  the  one  produced  by  the  freestream.  The  one  produced 
by  the  source  distribution  is  to  cancel  this  velocity  component  and  is  equal 
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to  equation  (49)  with  the  sign  reversed.  It  follows  that  the  source  strength 
can  be  calculated  as 

3zt 

0  =  cos  a 

The  induced  velocity  produced  by  equation  (46)  is  calculated  through  a 
source  lattice  method  which  is  derived  in  Appendix  A. 

Experimental  measurements  show  that  on  wings  with  asymmetric  thickness 
distribution,  the  lower  surface  has  virtually  no  effect  on  the  vortex 
development  (ref.  21).  Therefore,  to  account  for  the  flow  turning  angle  at 
the  leading  edge,  the  angle  of  attack  will  be  reduced  by  the  upper  surface 
nose  angle  6  (ref.  22).  In  other  words,  in  equation  (45)  a  is  to  be 

IN  S 

replaced  by  (a  +  SM).  The  nose  angle  for  a  symmetric  wedge  configuration  is 

equal  to  half  of  the  included  angle.  For  a  conventional  airfoil  section,  the 

slope  calculated  at  5%  local  chord  position  near  the  leading  edge  was  used  in 

reference  22  as  the  nose  angle  and  the  same  value  has  been  used  in  the  present 

method  as  well  (Fig.  6). 

2.3.4.  Effect  of  Boundary  Layer  Separation 

By  boundary  layer  separation  it  is  meant  that  after  separation  there 
is  a  thick  wake  without  a  strong  "organized"  vortex  flow,  as  usually  present 
on  a  wing  of  moderate  to  high  aspect  ratio  and  small  leading-edge  sweep  angles 
at  high  angles  of  attack.  To  calculate  the  associated  complex  flow  field 
would  require  solving  the  Navier-Stokes  equations  with  appropriate  turbulent 
modeling.  Since  this  is  not  feasible  at  the  present  time  for  preliminary 
design,  a  method  based  on  utilizing  nonlinear  section  data  was  developed  (ref. 
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10).  This  method  was  based  on  the  concept  of  matching  the  near  field  lift  (or 
vorticity)  with  the  far  field  value  provided  by  the  lifting  surface  theory. 

Based  on  this  concept,  it  is  assumed  that  the  effect  of  flow 
separation  is  to  reduce  the  local  angle  of  attack  by  Aa  at  a  spanwise 
station  (ref.  10).  Therefore,  the  effective  angle  of  attack  at  any  spanwise 
section  becomes 

a  =  a  -  a.  -  a  -  Act  (51) 

e  n  i  o 

where  an  is  the  geometric  angle  normal  to  the  section  which  may  have 

dihedral,  a.  is  the  induced  angle  of  attack,  a  is  the  angle  of  zero  lift, 
i  o 

and  Aa  represents  a  reduction  in  an  to  be  calculated  due  to  viscous 
effects.  It  follows  that 


Since  C£(3-d)  *s  computed  with  an  inviscid  theory,  its  value  is  usually 

larger  than  c  if  Aa  ■  0.  Therefore,  f  is  usually  less  than  1.0.  In 

i(2-D) 

this  case,  a  geometric  angle  of  attack  (a')  which  produces  the  reduced  lift 


can  be  found.  That  is, 


sin  a'  *  f  sin  a 

n 

or, 

a'  m  sin  ^(f  sin  a  )  (55) 

n 

If  follows  that  Aa  in  equation  (51)  becomes 

Aa  =  a  -  a'  (56) 

n 

The  solution  is  obtained  iteratively  as  follows: 

1.  Assume  Aa  *  0. 

2.  Find  a^  from  equation  (53). 

3.  Calculate  f  from  equation  (54). 

4.  Determine  Aa  from  equations  (55)  and  (56). 

5.  Use  Aa  to  reduce  a  in  the  3-D  boundary  condition  to  determine 
CA(3-D)* 

6.  Repeat  steps  2  through  5  until  the  successive  total  lift  coefficients 
differ  by  less  than  0.5%. 

For  a  wing  with  weak  or  moderate  vortex  flow,  such  as  that  in  a  strake- 
wing  configuration,  iterations  for  vortex  filaments  will  start  only  after  the 
iterations  for  nonlinear  section  data  are  completed. 

This  method  is  not  applicable  to  those  configurations  which  are  domin¬ 
ated  by  vortex  flow,  such  as  the  F-106B  configuration. 
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2.3.5.  Effect  of  Vortex  Breakdown 


The  effect  of  vortex  breakdown  is  one  of  the  important  factors 

affecting  the  aerodynamics  of  high  performance  aircraft.  Since  there  is  no 

comprehensive  theory  to  predict  the  breakdown  location  and  the  residual  vortex 

strength  after  breakdown,  a  semi-empirical  formula  derived  from  a  least-square 

analysis  of  available  data  is  used  in  the  program. 

Lamar  (ref.  23)  observed  that  for  delta  wings,  the  angles  of  attack 

for  vortex  breakdown  at  the  trailing  edge  (a  )  is  related  to  the  leading- 

BD 

edge  suction  distribution  (c  ).  Let 

s 


c  *  c  c  /c  sin  a 
s  s 


(57) 


where  c  is  the  mean  aerodynamic  chord  (MAC).  Since  cs  is  proportional  to 
2  — 

sin  a>  it  follows  that  c  is  a  function  of  planform  and  Mach  numbers.  From 

s 

a  least-square  analysis  of  Wentz'  data  (ref.  24),  it  was  found  that 

BD 

expressed  in  degrees,  would  fall  on  a  single  curve  described  as  follows  (Fig. 
7a) 

°BD  *  9‘195  ”  23-734y£  +  60-810yJ  -  33.533y3  +  7.391y£  "  0.581^, 

if  y  <  2.5. 

I 

ofiD  -  38.0,  if  y^  >  2.5  (58) 

where  y  is  the  distance  from  apex  to  centroid  of  the  c  distribution  from 
i  s 

inboard  to  q  of  c  measured  along  the  leading  edge  and  referred  to 

s \ max j 

the  half  span.  This  is  illustrated  in  Figure  7a. 


-23- 


was 


The  progression  of  the  breakdown  point  on  delta  wings  at  a  >  aDn 

DU 

also  analyzed  based  on  Wentz'  data.  Although  there  was  considerable  scatter 
in  the  experimental  data  (Fig.  7b),  a  single  curve  based  on  a  least-square 
analysis  was  obtained  to  be 

Ax  -  0.457(Aa)  -  0.1615(Aa)2  +  0.0303(Aa)3  -  0.00237(Aa)4 

+  0.00009(Aa)5,  if  Act  <  8.0  deg.  (59) 

Ax  ■  0.5392  +  0.0226(Aa),  if  Aa  >  8.0  deg. 

where  Aa  *  a  -  a„„  and  Ax  is  the  nondimensional  x-distance,  referred  to 
BD 

the  root  chord,  from  the  trailing  edge  to  the  breakdown  point. 

It  is  known  that  the  vortex  strength  after  breakdown  is  reduced  but  not 
vanished.  In  the  present  method,  the  vortex  strength  outboard  of  a  station 
where  vortex  breakdown  occurs  is  multiplied  by  a  factor  k  to  represent  the 
residual  vortex  lift.  The  factor  k  is  determined  again  by  analyzing  Wentz' 
data  and  is  found  to  be  (ref.  25) 

k  -  0.131  +  0 . 384y^ ,  if  <  1.49 

k  -  0.951  -  0.208y^  +  0.028yJ,  if  1.49  <  y^  <  3.71 

k  -  0.5,  if  y^  >  3.72  (60) 

In  addition,  camber  (including  vortex  flaps)  will  affect  a  .  In  the 

present  code,  it  is  assumed  that  the  leading-edge  vortex  will  burst  only  when 
it  reaches  the  strength  corresponding  to  a  flat  planform.  The  strength  can  be 
represented  by  the  total  leading-edge  suction  integrated  to  cs(max)'  Cgc 
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and  Cg£  are  these  integrated  values  for  a  cambered  wing  and  the 
corresponding  flat  wing,  respectively,  then 


C  sin^  (a-a  ) 
sc  'o' 

C  ,  2 

sf  sin  a 


(61) 


from  which  a  can  be  found.  a  is  to  be  added  to 


“bd 


obtained  earlier 


o  o 

for  a  flat  wing  when  camber  is  present. 

Furthermore,  it  was  found  that  upper  surface  slopes  in  the  spanwise 
direction  due  to  thickness  distribution  will  affect  the  movement  of  the  burst 
point  (ref.  26).  To  model  this  effect,  it  is  assumed  that  the  local  angle  of 
attack  is  changed  by  an  amount  equal  to  the  spanwise  upper  surface  angle; 
i.e.,  tan  ^(dz/dy).  The  new  local  angle  of  attack  is  used  in  equation  (59)  to 
determine  the  vortex  burst  point  iteratively  because  dz/dy  is,  in  general, 
not  a  constant.  This  effect  was  found  to  be  quite  significant  for  a  highly 
cambered  wing,  such  as  the  F-106B  configuration. 

For  a  wing  with  a  sideslip  angle  f3,  the  angles  of  attack  for  vortex 
breakdown  on  the  windward  side  (afiR)  end  leeward  side  (aBL)  are 

calculated  as  follows: 

1.  Calculate  y^  from  the  symmetrical  suction  distribution  with  the 
leading  edge  sweep  angles  -  p  for  the  right  wing  and  +  p  for 
the  left  wing. 

2.  Assume  that  the  maximum  vortex  strength  before  breakdown  is  unchanged 

by  sideslip.  Because  of  increase  in  vortex  lift  on  the  right  wing  in 

a  positive  sideslip,  the  maximum  vortex  strength  would  be  reached  at  a 

lower  a  (i.e.,  a__)  than  a  on  the  right  wing.  Since  the  vortex 
BR  BD 

2 

strength  is  represented  by  cg  which  is  proportional  to  sin  a,  ctgR 
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and  a  can  be  obtained  by  solving  the  following  equation: 

BL 

2  2  2  2 
sin  aBD  *  sin  afiR  +  (±  2  sin  p  +  tan  sin  p)sin  (62) 

where'  is  for  at  . 

BL 

The  effect  of  yawing  on  vortex  breakdown  is  similar  to  that  due  to 

sideslip.  Since  the  equivalent  sideslip  (p^)  is  variable  along  the  leading 

edge,  the  average  of  p^  at  a  given  y  station  and  that  at  the  tip  is  used 

to  determine  breakdown  a.  In  the  program,  a_T  for  the  yawing  motion  is 

determined  by  a  linear  interpolation  between  a__  and  ol0  for  sideslip. 

In  a  rolling  motion,  the  local  angle  of  attack  on  the  right  wing  is 

increased  due  to  a  positive  roll.  It  seems  that  the  vortex  breakdown  angle 

would  be  decreased.  However,  rolling  tends  to  move  the  centroid  of  the  vortex 

lift  distribution  outboard  (i.e.,  to  increase  y^).  Therefore,  the  vortex 

breakdown  angle  is  increased.  In  the  present  code,  the  vortex  breakdown  angle 

for  the  rolling  motion  is  assumed  to  be  equal  to  a  for  the  symmetrical 

BD 

loading. 


2.3.6.  Forebody  Vortex  Separation 

It  is  well  known  that  a  fighter  forebody  may  generate  symmetric  or 
asymmetric  vortex  separation  at  high  angles  of  attack  and  zero  sideslip. 
Asymmetric  separation  will  produce  side  force,  and  hence,  yawing  moment  (ref. 
27).  Although  the  vortex  lattice  method  has  been  used  to  model  the  symmetri¬ 
cal  vortex  separation  on  a  body  of  revolution  (ref.  28),  there  is  no  three- 
dimensional  modeling  of  asymmetric  separation  at  zero  sideslip  available  for 
applications  in  preliminary  design.  However,  it  was  shown  in  reference  29  by 
a  slender  body  theory  with  two  discrete  vortices  that  asymmetric  vortex 
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separation  could  be  obtained  as  the  additional  solution  in  a  boundary  value 
problem.  This  additional  solution  may  generate  a  positive  or  negative  side 
force  with  equal  probability.  This  concept  is  now  included  in  the  present 
code  as  follows. 

Calculated  results  from  the  slender  body  theory  are  valid  in  the  near 
field  and  provide  the  strengths  and  initial  locations  of  two  forebody 
vortices.  The  locations  of  these  forebody  vortices  in  the  wing  downwash  field 
are  adjusted  at  the  same  time  as  the  wing  free  vortex  filaments  are  adjusted 
in  position.  The  slender  body  theory  will  be  applied  at  an  effective  angle  of 
attack  (ag)  at  several  stations  of  the  forebody  to  include  the  wing  upwash 
and  body  camber,  if  any. 

Since  a  forebody  is  typically  not  of  the  circular  cross  section,  a 
numerical  conformal  mapping  is  applied  at  computational  stations  to  map  each 
forebody  cross  section  into  a  circular  one  with  radius  ”a".  The  mapping 
function  used  is  (ref.  30) 


Z 


k 


"C  +I 

n=l 


A 


n 


n+1 

a 


n 

C 


iC 


(63) 


where  Z  represents  a  complex  variable  in  the  physical  plane  and  £  the 
circle  plane  (see  Appendix  B) . 

The  complex  potential  on  the  circle  plane  includes  images  of  vortices  and 
a  source  term  to  represent  the  effect  of  expansion  of  body  cross  sections  (see 
Fig.  8): 

a2  a2  F1  5  _  C1 

W(?)  -  -iV,  sin  ae(?  ~  ~  \  sin  g(C  +  — )  - j 

cL 

c  *  — 


+  vz> 


(64) 


C 
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where  5  and  5  are  complex  conjugates  of  5  and  5  ,  respectively.  The 
*  ^  1  z 

complex  velocity  is  given  in  equation  (B.8)  of  Appendix  B.  The  boundary 

conditions  to  be  satisfied  are  (ref.  29): 

1.  At  separation  points,  s^  and  S2  (see  Fig.  8),  tangential  velocities 
corresponding  to  the  unknown  vortex  strengths  are  imposed. 

2.  The  vortices  are  force-free. 

The  resulting  nonlinear  algebraic  equations  are  solved  simultaneously  for  r  ^ , 
xl»  yi»  r2,  x2,  y2»  where  (x,y)  are  the  vortex  core  coordinates  in  the 

physical  plane.  Typically,  two  pairs  of  realis  ic  solutions — one  symmetrical 
and  one  asymmetrical — can  be  obtained.  The  asymmetrical  solution  produces  a 
side-force  and  hence,  a  yawing  moment  at  zero  sideslip.  In  lateral- 
directional  motion,  or  with  asymmetrical  separation  points,  two  branches  of 
solutions,  both  asymmetrical  can  be  obtained.  Since  it  is  not  known  when  the 
second  branch  (which  is  more  asymmetrical)  may  appear  under  a  given  flight 
condition,  the  second  branch  of  solutions  may  be  used  to  calculate  the  maximum 
possible  side  force  and  yawing  moment  for  design  purpose. 

Detailed  expressions  for  the  boundary  value  problem  in  slender  body 
theory  can  be  found  in  Appendix  B  and  reference  29. 


2.4.  Calculation  of  Aerodynamic  Characteristics 

The  pressure  distribution  on  lifting  surfaces  is  calculated  by  applying 

Kutta-Joukowski's  theorem.  There  are  two  types  of  contribution  to  the  lifting 

pressure  coefficients.  One  is  from  the  bounded  vortex  elements  and  the  other 

from  the  chordwise  vortex  elements.  aC  from  the  bounded  elements  can  be 

P 

written  as 


(AC  )  =  2(u  -  v  tan  $)y 

P  a 


(65) 


-28- 


where  u,v  are  the  total  x-  and  y-velocity  components  and  4;  the  element 
sweep  angle.  The  lifting  pressure  coefficient  from  a  chordwise  vortex  element 
has  the  following  form: 

(AC)  =  2Tv  (66) 
P  T 


Note  that  the  sidewash  (v)  is  mainly  produced  by  the  free  vortex  filaments. 

The  calculated  AC  is  interpolated,  if  necessary,  to  obtain  AC  at 

P  P 

integration  stations  and  may  be  printed  out  at  pre-selected  locations.  For 
more  details,  reference  11  should  be  consulted. 

The  calculated  pressure  force  is  assumed  to  be  acting  normal  to  the 
camber  surface .  This  pressure  force  is  then  resolved  and  integrated  to  obtain 
the  following  sectional  characteristics: 


,  Xte  3z  dz  2  5z  2  1/2 

c*  ”  c  /  iCpt  sin  a  +  cos  aVil  +  (q~)  +  iof-)  ]  d*  (67) 

XAe 

,  xte  az  az  2  az  2  1/2 

cd  =  c  /  ACpt'  aT-  cos  a  +  sin  +  ^aF-)  +  (^r)  1  dx  (68) 

XJLe 

,  Xte  az  dz  2  3z  2  1/2 

ACp(Ax  +  A*  +  (sr)  +  (#)  1  <69> 

C  XJle 


In  an  attached  flow,  c^  is  reduced  with  the  leading-edge  thrust.  The  total 
force  coefficients  are  calculated  by  spanwise  integration  of  sectional  force 
coefficients  as 


b/2 


X  '  S  [  Ct  Cdy 


2  b/2 

CD  "  S  £  Cd  Cdy 


(70) 


(71) 
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(72) 


„  2  (  2, 

C  =  —  I  c  c  dy 

m  —  J  m 


To  calculate  the  pressue  coefficient  on  the  fuselage,  the  tangential 
velocity  component  on  the  body  surface  is  needed.  For  this  purpose,  the  unit 
tangent  vector  on  the  surface  will  be  determined  first.  Note  that  in  a 
crossflow  plane,  the  unit  normal  vector,  n,  is  (see  Fig.  3  and  equation  (11)) 

♦  .  *  1  6R 

p  t  e  —  — 

->  r  9'  R  50 
/  1  +  ^R  50  J 1 

It  follows  that  the  unit  tangent  vector,  t,  is 

>  ,1  5R>,  •+■ 

t.ikael_LV  (73) 

A  +  fIMi2 

7  1  +  ^R  50^ 

With  the  total  velocity  given  by  equation  (13),  the  total  tangential  velocity 
in  a  crossflow  plane  is  then  given  by 


v  *  V  •  t 


For  example,  in  a  symmetrical  flight  condition, 


vt  *  [<vr  +  V.  sin  «  cos  0)  £  -  v0  -  V,  sin  a  sin  @]//l  +  [j  ||)2 


The  pressure  coefficient  is  then  (ref.  9) 


1  -  (1+u)2  +  M2  u2  +  (cos  a  |^)2  +  v2 

OO  v  5x '  t 
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The  sectional  normal  and  axial  force  coefficients  are  obtained  by 
decomposing  the  pressure  force  into  the  appropriate  directions  as  follows 


1  n 
-T—  /  C 


ref  —it 


p(f)  ]W|  * 


'C*  1  vf  ^  J  ,1  dR-,2 

^  **nj  =  "  Rref  [n  Cp(f)  j  VF  |  k  ^  +  (r  60 )  RdG 


1  % 
— J  c 


ref  —it 


P(f)  iwy 


•  tA*  (±  ||)2  Rde 


Using  equation  (11)  for  VF  and  equation  (12b)  for  k,  it  is  obtained  that 


ref  —it 


0  .  1  5R 

cos  8  +  sin  e  I  ^ 

A  +  (iH)2MS2 


i  +  (I  H)2  Rd* 


_  4 _ fC  _ _ ax 

Rref  -x  P(f)  /  /I  dR^2  .  r5R^2 

/  1  +  In  an  J  + 


1  +  H)2  RdG 


The  side  force  coefficient  can  be  calculated  in  a  similar  manner: 


cy  -  -  5—  /’  cP(£)  ifr  J  71  +  (5 1)2  Rde 


ref  —it 


-  -  ^  /’  cPCO  r  6  ;  IV  1 !  i  -  (30 

/L  AR59)  + 

By  integrating  sectional  coefficients  along  the  fuselage,  the  total  force 
coefficient  can  be  obtained.  For  example,  the  normal  force  coefficient  is 
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1  ft 

C „/,v  ■  r  f  R  ^  c  , . .  dx 
N(f)  S  J  ref  n(f ) 

su 

The  pitching  moment  coefficient  is 

i  *£t 

C»(£)  *  '  "=  I  Rre£  Cn(f  >  X<1X 

SC  Xfi 

3.  NUMERICAL  RESULTS  AND  DISCUSSION 

For  the  purpose  of  identification,  the  present  code  will  be  called  the 
VORSTAB-II  code  in  the  following.  In  using  the  code,  whenever  the  airfoil 
sectional  data  are  needed,  they  are  generated  with  a  modified  Eppler's  code 
(ref.  33).  In  the  following,  calculated  results  for  F/A-18,  F-5,  F-16XL,  and 
F-106B  configurations  will  be  presented.  All  calculations  are  made  in  eight 
iterations  to  adjust  the  free  vortex  and  wake  positions,  except  for  the  F16XL 
where  10  iterations  are  needed.  The  vortex-breakdown  effect  is  applied  only 
at  the  last  Iteration.  Results  for  a  less  than  10  deg.  are  all  obtained  with 
the  method  of  suction  analogy.  The  main  emphasis  is  on  lateral-directional 
aerodynamics. 

An  F/A-18  Configuration 

This  is  a  clean  configuration  without  deflection  of  the  leading-edge  flap 
as  shown  in  Fig.  9.  On  the  strake,  it  is  assumed  that  the  aerodynamics  are 
dominated  by  vortex  flow,  so  that  no  section  data  are  used.  On  the  other 
hand,  both  vortex  flows  and  viscous  separation  are  assumed  to  exist  on  the 
wing  at  high  angles  of  attack. 


(82) 


(83) 
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Longitudinal  characteristics  are  compared  with  data  in  Figure  10.  The 
results  indicate  that  C^,  Cp,  and  Cm  can  be  reasonably  predicted.  However, 
for  a  greater  than  20  deg.,  there  is  some  discrepancy  in  the  predicted 
results.  This  is  due  to  uncertainty  in  integrating  the  fuselage  loading. 
Since  the  inviscid  fuselage  loading  is  significant,  total  lift  would  be  over¬ 
predicted  if  all  fuselage  loading  is  included  in  the  calculation.  To  solve 
this  problem,  the  idea  used  in  the  USAF  Datcom  (ref.  35)  is  to  exclude  the 
region  of  separated  flow  from  lift  calculation.  However,  it  is  not  known 
where  the  fuselage  separation  should  occur.  In  the  present  coue,  it  is 
arbitrarily  set  at  a  location  where  the  sectional  lift  coefficient  is  1.5. 

Lateral-directional  coefficients  (C  ,  C  ,  C  )  for  this  configuration 

yB  £B  n8 

are  presented  in  Figure  11.  These  characteristics  are  significantly 

influenced  by  free  vortex  positions,  bursting  locations  and  forebody  vortices. 

The  vortex  system  is  illustrated  in  Figure  12.  In  a  positive  sideslip,  the 

right  leading-edge  vortices  are  closer  to  the  right  vertical  tail  than  the 

left  ones  to  the  left  vertical  tail.  At  an  angle  of  attack  of  25  deg.,  the 

strake  vortices  begin  to  burst  at  a  position  near  the  vertical  tail  (ref.  36). 

At  the  same  time,  the  forebody  vortices,  although  still  weak,  tend  to  push  the 

right  strake  vortices  towards  the  right  vertical  tail.  It  is  a  combination  of 

these  factors  that  generate  more  downloads  on  the  right  tail  than  on  the  left 

one  to  produce  a  negative  C  .  As  the  angle  of  attack  is  further  increased, 

n8 

the  vortices  will  move  upwards  to  reduce  the  interference  effect.  Overall, 
the  calculated  results  agree  well  with  data.  It  should  be  noted  that  the 
method  of  suction  analogy  is  not  capable  of  predicing  this  interference  effect 
(ref.  9).  The  empirical  prediction  of  vortex  breakdown  location  at  a  *  25 
deg.  has  been  shown  to  be  in  good  agreement  with  water-tunnel  results  (ref. 
36). 
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An  F-5  Configuration 


The  F-5  configuration  is  particularly  notable  for  its  long  forebocy  (fig. 

13).  The  whole  wing  is  assumed  to  be  subject  to  the  effects  of  viscous 

separation  and  vortex  flow.  The  results  for  lateral-directional  coefficients 

are  presented  in  Figure  14.  The  calculated  results  indicate  the  loss  of 

directional  stability  at  a  slightly  lower  angle  of  attack  than  the  data 

shown.  The  more  noteworthy  is  the  effect  of  forebody  vortices.  The 

asymmetric  forebody  vortices  at  a  positive  sideslip  are  known  to  produce  a 

positive  sideforce  to  contribute  to  a  position  directional  stability, 

i.e.,  +C  (ref.  38).  As  it  is  shown  in  Figure  14,  this  effect  of  forebody 
n8 

vortices  becomes  significant  at  a's  above  25  deg.  to  produce  a  positive  C 

n 


Note  that  typically  a  vertical  tail 
high  a's.  The  increase  in  dihedral 


will  not  produce  a  positive  C  at 

nB 

effect  (C  )  at  higher  a's  is  mainly 
*8 


produced  by  an  assumption  in  the  code  that  the  loading  near  the  left  wing  tip 
is  reduced  in  proportion  to  the  sideslip  angle  because  it  behaves  more  like 
the  wing  trailing  edge  as  the  sideslip  angle  is  increased.  In  the  meantime, 
the  right  wing  tip  will  have  its  loading  increased  to  simulate  the  leading- 
edge  effect. 


An  F-16X.L  Configuration 

Data  for  this  configuration  are  taken  from  ref.  39.  It  does  not  have 
camber  and  has  an  inboard  sweep  angle  of  70  deg.  and  an  outboard  sweep  of  50 
deg.  (Fig.  15).  Because  of  the  expected  strong  vortex  flow,  it  is  assumed 
that  it  is  dominated  by  vortex  flow  and  that  no  airfoil  section  data  should  be 
used. 
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The  lateral-directional  g-derivatives  are  presented  in  Figure  16.  It  is 

seen  that  C  is  overpredicted  at  a  of  35  deg.  The  most  difficult  task,  in 

£S 

modeling  this  configuration  is  again  calculating  the  fuselage  loading. 
Because  the  strong  vortex  flow  tends  to  produce  a  large  loading  on  the 
fuselage,  a  separated-flow  model  for  the  fuselage  would  be  needed  for  a  better 
prediction  of  C  .  Although  the  predicted  C  variation  is  reasonable,  the 

ye  ne 

calculated  C  is  more  erratic.  In  the  present  calculation  of  lateral- 

y8 

directional  characteristics,  the  fuselage  loading  is  integrated  only  from  the 
nose  to  a  station  at  which  the  sectional  lift  coefficient  is  less  than  1.0  or 
the  rate  of  change  of  sectional  lift  coefficients  is  less  than  1.0. 


An  F-106B  Configuration 

The  wing  has  a  conical  camber  and  a  leading-edge  sweep  angle  of  60  deg. 

(fig.  17).  It  is  assumed  to  be  dominated  by  vortex  flow  so  that  nonlinear 

sectional  data  are  not  used.  The  forebody  is  modeled  with  a  distribution  of 

vertical  elliptical  cross  sections.  The  forebody  is  transitioned  into 

circular  cross  sections  in  the  afterbody. 

A  typical  distribution  of  vortex  filaments  for  this  configuration  in  a 

positive  sideslip  is  illustrated  in  Figure  18.  Again,  at  high  a's, 

uncertainty  in  fuselage  modeling  in  both  geometric  shape  and  loading 

integration  causes  C  to  become  too  negative  at  a  *  30  deg.  (fig.  19).  At 

yB 

a  =  35  deg.,  the  predicted  C  is  not  negative  enough.  Again,  this  may  be 

ne 

related  to  difficulty  in  fuselage  modeling.  As  shown  in  this  figure,  C 

e 

appears  to  be  well  predicted. 

Steady  roll-rate  derivatives  for  this  configuration  are  presented  in 
Figures  20  and  21.  The  present  results  are  quite  similar  to  those  given  by 
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the  method  of  suction  analogy  (ref.  9)  Good  agreement  with  data  can  be 
obtained  before  vortex  breakdown.  However,  large  discrepancy  between  data  and 
theoretical  results  exist  at  high  a*  This  is  because  the  experiment  data 
were  obtained  in  forced-oscillation  tests.  Therefore,  unsteady  aerodynamic 
phenomena,  such  as  vortex  lag  and  unsteady  vortex  bursting,  would  become 
important.  Specifically,  C  must  be  included  in  the  prediction.  A  steady 

£e 

flwo  model,  such  as  the  present  one,  cannot  account  for  such  phenomena. 


4.  CONCLUSIONS 

A  steady-flow  aerodynamic  method  was  developed  to  calculate  lateral- 
directional  aerodynamics  at  high  angles  of  attack.  In  the  method,  the  edge- 
separated  vortex  flow  was  represented  by  free  vortex  filaments  and  visous 
separation  effect  was  accounted  for  by  using  nonlinear  section  data.  Forebody 
vortex  separation  was  calculated  with  a  slender  body  theory  at  an  effective 
angle  of  attack  which  was  given  by  the  lifting  surface  theory.  Calculated 
results  for  four  configurations  were  compared  with  available  wind-tunnel  data. 
Reasonable  agreement  could  be  obtained  for  sideslip  derivatives.  Discrepancy 
in  C  was  attributable  to  inaccurate  fuselage  modeling.  It  was  also  clear 

ye 

that  dynamic  stability  coefficients  at  high  angles  of  attack  could  not  be 
accurately  calculated  with  a  steady  flow  theory. 
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APPENDIX  A 


A  SURFACE  SOURCE  LATTICE  METHOD 


Assume  that  the  continuous  source  distribution  is  reduced  to  a  finite  sum 
of  discrete  source  element  contributions .  These  source  elements  occupy  the 
same  locations  as  the  vortex  lattices.  A  source  element  is  defined  by  its 
inboard  endpoint  coordinates  (x^,y^,z^)  and  outboard  endpoint  coordinates 
(x2,y2>z2)*  A',  integration  parameter  x  *  representing  fraction  of  the 

element,  can  be  defined  such  that 

x'  -  x  =  x  -  x  +  t  (x2  -  x1 ) 

y'  ~  y  *  yL  -  y  +  t(y2  -  yL)  (a.i) 

z'  -  Z  -  Z  -  2  +  x(z2  -  Zx) 

Therefore,  the  denominator  in  the  integrand  of  equation  (46)  can  be  written  as 

(x'-x)2  +  s2(y'-y)2  +  g2(z'-z)2  =  ^[(x^)2  +  g2(y2-y1)2 

+  62(z2_z1)2]  +  2t[(x1-x)(x2-x1)  +  62(y1-y)(y2~y1) 

+  g2(z1-z)(z2-z1)]  +  (xL-x)2  +  g2(y1~y)2  +  62(z1~z)2 
-t2  A  +  tB  +  C 

where 
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(A. 2a) 


(x2~xi )2  +  82(y2_y1)2  +  82(z2_zi)2 

2[  (x1~x)(x2-x1)  +  82(y1-y)(y2~y1)  +  e2(z1-z)(z2-z1>] 

(x^x)2  +  g2(y1-y)2  +  g2(z1-z)2 


(A. 2b) 
(A. 2c) 


It  follows  that  equation  (46)  can  be  reduced  to  a  sum  of  quantities  of  the 
following  form: 


oAA  r 
Un  J 


0  /at2  +  Bt  +  c 


2— - | £n[2A  +  B  +  2/ A  /A  +  B  +  C 


-  Jtn[B  +  2/ AC]} 


7  2  2  2 

A£  =  /(x2-x1)  +  (y2-y1)  +  (z2~z^) 


(A. 3) 


By  differentiation,  the  following  induced  velocity  components  can  be  obtained: 

-  xo  _  xi 

x„-x,  +  2/A  — - - —  .  „  7a  , 

...  21  /A  .  p  ,  r  x,  -  x,  +  2  / -z  (x  -x) 

3A$  m  ctAA  r _ / A  +  B  +  C  _  __2 _ 1 _ /  C  1  |  (A. 4) 
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(A. 5) 
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APPENDIX  B 


CALCULATION  OF  FORE BODY  FLOW  FIELD  IN  THE  CROSSFLOW  PLANE 


By  slender  body  theory,  aerodynamic  characteristics  of  a  forebody  with 
vortex  separation  can  be  calculated  in  the  crossflow  plane.  The  body  cross 
section  is  first  mapped  conformally  to  a  circle  of  radius  "a."  The  complex 
potential  on  the  circle  plane  can  be  easily  constructed.  The  total  force 
coefficients  can  be  obtained  by  contour  integration.  These  are  summarized  in 
the  following. 


B.l  Conformal  Mapping 

In  reference  30,  a  formula  for  mapping  an  arbitrary  polygon  symmetrical 
about  the  real  axis  to  a  circle  of  radius  "a"  is  given: 


k  A  a 

Z1  '  +  l,  — 

n-1  C]_ 


(B.l) 


where  C  is  a  constant  used  to  adjust  the  location  of  the  center.  Since  the 
cross  section  of  a  forebody  is  usually  symmetrical  about  the  vertical  axis,  a 
rotation  by  tt/2  is  needed  (see  Fig.  B.l): 


(B.2) 


It  follows  that  equation  (B.l)  becomes 


k  A  a 

?  +  I  — 

n-1  5" 


(B.3) 


The  derivative  dZ/d^  is  needed  in  mapping  of  flow  nets: 


k  n+1  , 

,  r  .  a  1 

1  -  )  n  A  — r: - tt 

L ,  n  n+1  ,  n+1 
n=l  5  1 


k  nA 


1  -  V  —  - 

1  L  *n+l  n+1 

n=l  5  i 


(B.4) 


where 


(B.5) 


For  an  elliptic  cross  section,  with  a^  and  &2  as  major  and  minor  axes, 
an  analytical  form  of  the  mapping  function  is 


Z  -  5  + 


2  2 
al  ~  a2  1 


(B.6) 


Comparing  with  equation  (B.3),  it  can  be  shown  that 


2  2 
al  S2 


*1  +  a2 


A  =0,  n=2,3,»»» 

n 


B.2  Boundary  Conditions  with  Vortex  Separation 


The  complex  potential  can  be  written  as 


2  2  r.  ? 

-iv.  sin  «(c  -  7")  '  v«  sin  &U  +  f")  +  27T  in  - 


'l  4  ~  ?2 

i7T*n - 2  +  WS(Z) 

a 

c  -  — 

^  2 
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(B.7) 


The  complex  velocity  is  given  by 


1  dW 


V  sin  a  d^ 


C  ?  C  “?]_ 


ihrl 


+  - V"]  +  V~sinT 


5  “5,  =* 


dW 

C 


4c 

5  " 

«1 


(B  .8) 


where 


1  2u  V  a  sin  a 

CO 


(B.9a) 


Yo  “ 


2  2ir  V  a  sin  a 


(B.9b) 


In  equation  (B.8),  the  source  term  W  is  to  simulate  the  flow  field  due 

s 

to  change  in  body  cross  sections  in  the  x-direction.  This  source  flow  must 
satisfy  the  following  condition: 


V 

CO 


cos  a 


dR 

dx 


(B.10) 


In  a  crossflow  plane,  the  source  flow  can  be  simulated  by 


$  =  C.  (x,0  )ju  r 

s  l 


(B. 11) 


where  is  regarded  as  a  function  of  0  because  of  R(x,0).  Therefore, 


CL(x,0) 


cos  a)RS 


(B.12) 


The  induced  velocity  due  to  the  source  flow  is  computed  from 
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(B. 13) 


3$c  C  (x,e)  (V  cosa)R^ 
si  oo  ax 

'  s  -  ss  ■ 

3  r  r  r 


Note  that  on  the  body  surface,  r  =  R.  in  addition,  for  a  flat  surface,  such 
as  a  delta  wing,  the  limiting  form  of  equation  (B.13)  at  the  leading  edge 
would  be 


34> 

— -  =  V 
3y  » 


cos  a 


d£ 

dx 


V  cos  a  tan  6 

00  X 


(B. 14) 


where  6x  is  the  half  apex  angle  of  the  delta  wing.  The  velocity  given  by 
equation  (B.14)  is  used  in  setting  up  the  conditions  at  the  separation  points 
which  are  sharp-edged,  such  as  on  a  delta  wing  or  a  chine  cross  section. 

The  separation  locations  are  calculated  with  Stratford's  criteria  in  the 
same  way  as  used  in  reference  31,  based  on  the  velocity  distribution  computed 
in  the  three-dimensional  flow  model,  including  wing  effect.  At  the  separation 
points  and  the  vortex  locations,  the  following  conditions  must  be  satisfied. 

1 .  Conditions  at  Separation  Points 

On  a  smooth  surface,  the  mean  tangential  velocity  (Vf  )  at  the 
separation  points  in  the  physical  plane  is  equal  to  the  convection  velocity  of 
vortex  in  the  shear  layer  after  boundary  layer  separation  (ref.  29).  It  is 
taken  to  be  (referring  to  eqns.  (3)  and  (7a)  of  ref.  29) 
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In  equation  (B.15),  K  is  taken  to  be  0.6  to  account  empirically  for  secondary 

separation  which  reduces  the  strength  of  the  primary  vortex  when  only  the 

latter  is  used  in  the  model,  f  is  taken  to  be  0.61  to  account  for  the  vortex 

convection  speed  in  a  shear  layer  (ref.  29).  This  tangential  velocity  (Vr  ) 

cm 

should  be  equated  to  the  total  tangential  velocity  calculated  in  the  flow 
model : 


k  dW, 
dZJ 


or 
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V  sin  a 


=  -  Im{  i 


10. 


V  sin 


dW 


d£_l 

dZJ 


at  Z  =  Z 


(B. 16) 


sep 


where  represents  a  tangent  angle  and  is  defined  in  Figure  B.2,  and  Im 
stands  for  the  imaginary  part.  Using  equation  (B.8),  equation  (B.16)  can  be 
rewritten  as 
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(B.17) 


* 

where  e  is  0,  in  the  circle  plane. 

When  Zaan  is  a  sharp  edge  so  that  the  mapping  is  singular, 

(dW/dj;  )(d$/dZ)  must  be  first  reduced  by  L'Hospital's  rule.  The  resulting 
velocity  will  consist  of  a  tangential  and  a  normal  component.  The  tangential 
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component  is  used  in  equation  (B.16),  but  the  normal  component  must  be  set  to 

zero.  This  introduces  one  extra  equation.  In  this  case,  f//K  in  V  must 

m 

be  regarded  as  unknown  to  result  in  a  determinate  set  of  equations. 


2 .  Force-Free  Conditions 

The  separated  vortex  should  be  force-free,  or  based  on  Kutta-Joukowski's 
theorem. 
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where  ry  an d  Zy  may  be  taken  to  be  r^,  or  and 

dt  by  dx/V  cos  a  and  considering  F  =  0,  it  is  obtained 

oo  V 

side  vortex  that 
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(B. 19) 


Using  the  complex  potential  in  equation  (B.8),  equation  (B.19)  can  be  reduced 
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(B.20) 


The  same  consideration  applies  to  the  left-side  vortex. 

For  a  body  with  a  smooth  surface,  there  are  two  equations  representing  the 
separation  condition  and  four  equations  (one  being  the  real  and  the  other 
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imaginary  part  of  equation  (B.20))  from  the  force-free  conditions.  These  six 
equations  are  to  be  solved  for  y^ ,  z^  ,  y^ ,  and  z^.  These  equations 
are  highly  nonlinear  and  they  are  solved  by  a  modified  Newton  iterative  method 
(ref.  29).  Typically,  a  symmetric  vortex  solution  can  be  obtained  from  any 
initial  input  data.  To  search  for  any  additional  solutions,  a  method  of 
function  deflation  to  remove  the  symmetric  vortex  solution  is  used  to  avoid 
convergence  to  the  same  symmetric  solution  again  (ref.  29).  The  additional 
solution  show  asymmetric  vortex  locations  and  strengths. 

B.3  Force  Coefficients  by  Momentum  Theory 

With  the  vortex  locations  and  strengths  determined,  the  complex  potential 
may  be  used  to  calculate  the  velocity  on  the  body  surface,  and  hence  the 
pressure  distribution.  The  latter  may  be  integrated  to  obtain  the  force 
components.  Since  the  force  components  are  of  primary  interest  in  the  present 
applications,  an  alternate  method  to  calculate  them  is  to  apply  the  momentum 
equation.  Based  on  Ward  (ref.  32),  the  lateral  force  based  on  slender  body 
theory  is  given  by 

- 1 -  -  1-21  6  2~  dZ  -  21  S„(sin  e.  -  1  sin  0)1  (B.21) 

(  V2  )pv„2  Sre£  sref  v-  B 

where  Sfi  is  the  base  area.  The  contour  integral  in  equation  (B.21)  can  be 
enlarged  so  that  both  W(^)  and  dZ/d^  can  be  expanded  for  large  £.  Using 
equations  (B.4)  and  (B.7),  it  can  be  shown  for  large  5  that 
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The  contour  integral  can  therefore  be  calculated  by  residue  theorem  as 
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Using  equation  (B.21),  the  force  coefficients  can  be  obtained: 
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the  force  coefficients  become 
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Figure  B.l.  Sequence  of  Transformations. 


Figure  B.2  Velocity  Components  at  a  Separation  Point  S^- 
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Figure  1.  Illustration  of  Flow  Features  Predicted 
in  the  Present  Method 
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Figure  2.  Difinition  of  Airplane  Coordinate  Systems 
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Figure  3.  Definition  of  Variables  Defining  a  Nonplanar 
Wing  Surface  (Ref.  13) 
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O  Control  Point 


A  Typical  Vortex  Element  of  Leading-Edge 
Vortex  System. 


Figure  6.  Effective  Angle  of  Attack  for  Thick 


(a)  Typical  Leading-Edge  Suction  Distributions 
of  Delta  Wings  with  Sharp  Leading  Edges 


Figure  7.  An  Empirical  Model  for  Vortex 
Breakdown  Effect  (Ref.  25) 


a  Della  Wings 
o  Cropped  Delta  Wings 
e  Cropped  Arrow  Wings 
—  Least  Square  Average 


Figure  7.  Concluded. 


Figure  9.  A  Geometric  Model  for  an  F/A-18 
Configuration 
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DA1A  (Ref .  34) 


0.2,  IL=  0.6x10 


Sketch  of  an  F-5  Configuration  (Ref.  37) 
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Figure  14.  Sideslip  Derivatives  of  an  F-5  Configuration 

z: 

at  M  =  0.1  and  0.56x10  Based  on  Body  Axes. 
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Three-View  Sketch  of  an  F-16XL  Model  (Ref.  39) 


Figure  17.  Three-View  Sketch  of  an  F-106B  Model  (Ref. 
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Figure  18.  Vortex  System  for  an  F-106B  Configuration  at 
an  Angle  of  Attack  of  25  deg.  and  M  =  0.1. 
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